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Abstract 

Three core-softened families of potentials are checked for the presence of density and diffusion 
anomalies. These potentials exhibit a repulsive core with a softening region and at larger distances 
an attractive well. We found that the region in the pressure-temperature phase diagram in which 
the anomalies are present increases if the slope between the core-softened scale and the attractive 
part of the potential decreases. The anomalous region also increases if the range of the core- 
softened or of the attractive part of the potential decreases. We also show that the presence of the 
density anomaly is consistent with the non monotonic changes of the radial distribution function at 
each one of the two scales when temperature and density are varied. Then, using this anomalous 
behavior of the structure we show that the pressures and the temperatures in which the radial 
distribution functions of the two length scales are equal are identified with the Widom line. 

PACS numbers: 64.70.Pf, 82.70.Dd, 83.10.Rs, 61.20.Ja 
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I. INTRODUCTION 



Core-softened {CS) potentials have been attracting attention due to their connections 
with the anomalous behavior of liquid systems including water. These potentials, U(r), 
exhibit a repulsive core with a softening region limited by r% < r < r 2 where d(rf)/dr > 
with / = —dll/dr—. Despite their simplicity, these models originate from the desire of 
constructing a simple two-body isotropic potential capable of describing the complicated 
features of systems interacting via anisotropic potentials ——. This procedure generates 
models that are analytically tractable and computationally less expensive than the atomistic 
models. Moreover, they are lead to conclusions that are more universal and are related to 
families of atomistic systems^— . 

One of the features that has been successfully described by many of these models is the 
density anomaly. For water the specific volume at ambient pressure starts to increase when 
cooled below T 4°C. The anomalous behavior of water was first suggested 300 years ago^ 
and was confirmed by a number of experiments^ 1 ^. Besides, between 0.1 MPa and 190 MPa 
water also exhibits an anomalous increase of compressibility— >2i and, at atmospheric pres- 
sure, an increase of isobaric heat capacity upon cooling^^. For the case of water the density 
anomaly is attributed to the presence of hydrogen bonds between neighbor molecules. As the 
temperature increases the bonds break and the density increases. However, other systems 
such as Te,— Ga, Bi,— S,— Gei 5 Te 85 ,— , silica,—^— silicon 4 ^ and BeF 2 ,— show the same 
density anomaly without presenting hydrogen bonds what suggests that the mechanism for 
the presence of density anomaly might be more universal. 

In compass with the presence of the density anomaly in water a few years ago it was 
suggested that there are two liquid phases, a low density liquid (LDL) and a high density 
liquid (HDL)^. The critical point ending this transition, found only in computer simulations 
is located at the supercooled region beyond the line of homogeneous nucleation and thus 
cannot be experimentally measured. Even with this limitation, this hypothesis has been 
supported by indirect experimental results^ 1 ^ 1 ^. The presence of two liquid phase and of 
second critical point is also observed in certain CS potentials^"— i22r^Z. 

Which are the conditions for a CS potential to exhibit density anomaly and two liquid 
phases? A definitive answer to this question is still missing. There are, however, a few 
clues. If a CS potential has discontinuous forces it presents two liquid phases but no density 
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anomaly 4 ^ is observed. However, once the CS potential is modified to have continuous forces, 
the anomalies appear—. 

Recently it has been proposed that a CS potential exhibits density anomaly if the two 
length scales identified with the softened region would be accessible^ 1 ^ 4 - 1 ^ what can be un- 
derstood as follows. The radial distribution, g(r) of a CS potential has peaks at g{r\) 
and g{r2) where 7*1 and r2 > r% are the two length scales of the CS potentiate. If 
dg(r) /dp\ r=ri dg(r) /dp\ r=r2 < 1 then the system would have density anomaly. 

In this paper we test if this link between the behavior of the structure (radial distribu- 
tion function) and the thermodynamic anomalies holds for a number of two length scales 
potentials. We study the pressure temperature phase diagram of a two Fermi model 4 ^. The 
advantage of this model is that by changing few parameters is possible to vary the distance 
and the difference in energy between the two length scales without introducing extra scales. 
Moreover the length scales are well defined. 

Hence, having identified a connection between the density anomaly and the behavior of 
the structure, we also test if the radial distribution function is also related to the presence 
of two liquids predicted for these CS potentials. We show that the pressures and temper- 
atures in which the radial distribution function associated with one scale equals the radial 
distribution function of the other scales is linked with peaks in the constant pressure specific 
heat, namely the Widom line. 

The remaining of this paper goes as follows. In Sec. [Ill the model is introduced and 
the simulation details are presented. In Sec. [TTT] the pressure-temperature phase diagram is 
presented together with the behavior of the radial distribution function with density and 
temperature. Conclusions are presented in sec. [IVl 

II. THE MODEL 

Our system consists of N identical particles interacting through a continuous pair poten- 
tial obtained by the addition of 3 different Fermi-Dirac distributions 49 , 



The resulting expression describes a family of pair interaction potentials discriminated by 
different choices of the parameters {e^, r oi , a iy a{\. Appropriated choices of the parameters 



3 




(1) 



3 



allow us to obtain potentials that go from a smooth two length scales potential to a sharp, 
almost discontinuous, square potential^^. 

In the Table [I] nine different sets of parameters are shown, organized in three families 
named, S, A, and R. As shown in Figure [IJ for each family a specific characteristic of pair 
interaction potential is tuned. Then it is possible to test the effect of changing the two 
length scales in the pressure temperature phase diagram. 

In the potentials S the slope between the two length scales is varied. Then it is possible 
to check if the slope between the two length scales controls the location in the pressure- 
temperature phase diagram of the density anomalous region as suggested by Yan et al.— . 

In the case of the potentials A, the attractive length becomes broader. Consequently 
using this potential we test if increasing the range of the attraction leads to a decrease in 
the critical pressure as proposed by Skibinsky et al.-. 

In the case of the potentials R, the repulsive length scale becomes broader. Therefore this 
family of potentials is appropriated to observe if the enlargement of the repulsive length scale 
leads to a decrease in the liquid-liquid critical pressure and to an increase in the liquid-liquid 
critical temperature as suggested by Skibinsky et alA In addition to verify the assumptions 
of Yan et al.— and of Skibinsky et alA related to criticality, these three families of potentials 
are the perfect scenario to check our hypothesis that the density anomaly region in the 
pressure-temperature phase diagram is delimited by properties of the radial distribution 
function at the two length scales. 




r r r 



FIG. 1. Interaction potential. 

The thermodynamic and dynamic behavior of the systems were obtained using NVT 
molecular dynamics using Nose-Hoover heat-bath with coupling parameter Q = 2. The 
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system is characterized by 500 particles in a cubic box with periodic boundary conditions, 
interacting with the intermolecular potential described above. 

Standard periodic boundary conditions together with predictor-corrector algorithm were 
used to integrate the equations of motion with a time step At* = 0.002 and potential cut off 
radius r* = 2.5. The initial configuration is set on solid or liquid state and, in both cases, the 
equilibrium state was reached after t* eq = 1000. From this time on the physical quantities 
were stored in intervals of At* R = 1 during t* R = 1000. The system is uncorrelated after 
t* d — 10, from the velocity auto-correlation function, and 50 decorrelated samples were used 
to get the average of the physical quantities. The thermodynamic stability of the system 
was checked analyzing the dependence of pressure on density, by the behavior of the energy 
and also by visual analysis of the final structure, searching for cavitation. 

In what follows we take £i, a = r o i/0.950 as fundamental units for energy and distance, 
respectively, and all physical quantities are expressed in reduced units, namely 

T * _ k B T 
£i 

P* = P^ 



(2) 



where T, P and D are respectively temperature, pressure and diffusion coefficient. The 



TABLE I. Parameters for potentials S, A and R in reduced units of e and a = r o i/0.950. 



Parameter 



^2 



'ol 



' o2 



'o3 



a 



1,2,3 



a. 



2,3 



1.000 1.000 1.000 
3.000 3.000 3.000 
1.000 1.027 1.023 
0.950 0.950 0.950 
1.400 1.400 1.400 
1.950 1.950 1.950 
0.025 0.040 0.055 
0.000 0.000 0.000 
1.000 1.000 1.000 



A 1 A 2 



A, 



1.000 1.000 1.000 

3.000 3.000 3.000 

1.023 1.023 1.023 

0.950 0.950 0.950 

1.400 1.400 1.400 

1.755 1.8525 1.950 

0.025 0.025 0.025 

0.000 0.000 0.000 

1.000 1.000 1.000 



Rl i?2 R3 



1.000 1.000 1.000 
3.000 3.000 3.000 
1.023 1.023 1.023 
0.950 1.050 1.150 
1.400 1.400 1.400 
1.950 1.950 1.950 
0.055 0.055 0.055 
0.000 0.000 0.000 
1.000 1.000 1.000 
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diffusion coefficient is obtained from the expression: 

D = lim (fc(* + t ) -Q( t o)] 2 )'o ( 3 ) 

*-+oo 6t 

where fj(t) are the coordinates of particle j at time t, and (■ • • } to denotes an average over 
all particles and over all to- 

The error associated with pressure and temperature are Ap* pa 0.005 and AT* pa 0.01. 



III. RESULTS 

A. Pressure- Temperature Phase Diagram 

Fig. [2] presents the pressure versus temperature phase diagram obtained for the three 
families: S, A and R. In all the nine cases the system exhibits at high temperatures a fluid 
phase, at intermediate temperatures and very low pressures a gas phase and at intermediate 
pressures a low density liquid phase (LDL) while at very high pressures a high density liquid 
phase (HDL). The coexistence line between the gas and the low density liquid phases (not 
shown) ends in a gas-LDL critical point illustrated as a filled circle in Fig. The LDL-HDL 
coexistence line (not shown) ends in a LDL-HDL critical point also shown as a filled circle. 
The two critical points are located in the pressure and temperature phase diagram by the 
point in which the isochores meet. The critical pressures and the critical temperatures values 
are confirmed by analyzing the slope of the pressure versus density at constant temperature 
phase diagram. The maximum of these curves identify the critical point. For the other 
state points the slope of the pressure versus density phase diagram is also used as a check 
of stability. 

The Table [Til and the Table HTT1 and the Fig [3] summarize the values of the first (liquid-gas) 
and second (liquid-liquid) critical points and their changes in the p — T phase diagram for 
the three families studied. 

Fig. [2] shows that in the family of potentials S the values of the pressure and temperature 
of the liquid-gas and the liquid-liquid critical points are not sensitive to the change of slope 
as predicted by Yan et al.— . For the A family, also illustrated in Fig. |2] indicates that the 
temperature of the liquid-gas critical point increases with the increase of the range of the 
attractive scale, while the temperature and the pressure of the liquid-liquid critical point 
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FIG. 2. Pressure-temperature phase diagrams for the S, A and R families of potentials. The gray 
lines are isochores, the solid lines are the temperature of maximum density (TMD), the dashed 
lines are the extrema of diffusion and the dot-dashed line are the extrema of the translational order 
parameter. The filled circles are the liquid-gas (at high temperatures) and the liquid-liquid (low 
temperatures) critical points. 



decrease. This result indicates that if the attractive scale increases the high density liquid 
requires less pressure to be formed while the gas phases exist for higher temperatures as 
predicted by Skibinstky et al.->^. In the case R, shown as well in the Fig. |2j the liquid- 
liquid critical pressure decreases with the increase of the range of the repulsive scale. This 



TABLE II. Liquid-gas critical point location for potentials S, A and R. 



Potential 


K p* c1 


Potential 




Potential 


K P * C1 


Si 


0.05 1.00 


At 


0.04 0.68 


Ri 


0.02 0.90 


S 2 


0.05 0.99 


A 2 


0.03 0.82 


R2 


0.04 0.88 


s 3 


0.04 0.88 


A 3 


0.05 0.90 


Rs 


0.04 0.88 



TABLE III. Liquid-liquid critical point location for potentials S, A and R. 



Potential 


T* v* 

± c 2 fez 


Potential 


C2 


Pc 2 


Potential 


T* v* 


Si 


1.31 0.54 


Ai 


2.75 


0.53 


Ri 


1.34 0.44 


s 2 


1.26 0.48 


A 2 


1.98 


0.50 


R2 


2.67 0.37 


S3 


1.34 0.44 


A 3 


1.33 


0.44 


R 3 


6.01 0.26 



result indicates that as the repulsive scale becomes broader, it requires less pressure for the 
high density liquid to be formed while the repulsive scale has almost no effect in the low 
density liquid-gas coexistence line as predicted also by Skibinstky et al. 8 ' 52 . A summary of 
the liquid-gas and liquid-liquid critical pressures and temperatures are shown on Fig. |3j 

1 •— 

« K , case 
■ K , case . 
Rj case 



0.3 0.6 0.9 0.3 0.6 0.9 0.3 0.6 0.9 

* * * 

T T T 

FIG. 3. Location of the critical points on pressure-temperature phase diagram for cases S, A and 
R. 

B. Density, Diffusion and Translational Anomalies 

Fig. [2] shows the temperature of maximum density (TMD) for all the nine studied cases as 
a solid thick lines. For all the potentials S, A and R the TMD lines are observed. The limits 
of the TMD in the pressure-temperature phase diagram are shown in the Table HVl where pi 
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TABLE IV. Limit values for density (p*), temperature (T*) and pressure (p*) of the thermo- 
dynamics anomalies on pressure-temperature diagram. Here the point pi represents the density, 
temperature and pressure of the point of the lowest pressure in the TMD line, p m represents the 
point of the highest temperature and ph represents the point of the highest pressure of the TMD 
line. 



case pi p m p h 


case pi p m Ph 


case pi p m Ph 


p* 0.34 0.35 0.38 
Si T* 0.70 0.73 0.60 
p* 0.70 0.73 1.14 


p* 0.32 0.34 0.39 
S" 2 T* 0.68 0.80 0.60 
p* 0.65 0.96 1.27 


p* 0.31 0.34 0.39 
S3 T* 0.60 0.82 0.61 
p* 0.61 0.92 1.29 


p* 0.42 0.43 0.45 
A x T* 0.61 0.71 0.59 
p* 2.15 2.29 2.51 


p* 0.37 0.39 0.42 
A 2 T* 0.51 0.79 0.69 
p* 1.33 1.56 1.82 


p* 0.31 0.34 0.39 
A 3 T* 0.60 0.82 0.61 
p* 0.61 0.92 1.29 


p* 0.31 0.34 0.39 
Rx T* 0.60 0.82 0.61 
p* 0.61 0.92 1.29 


p* 0.33 0.37 0.42 
R 2 T* 0.77 0.92 0.50 
p* 1.28 1.81 2.67 


p* 0.34 0.39 0.47 
R 3 T* 0.70 0.92 0.32 
p* 2.04 3.65 5.76 



represents the values of (p*,T*,p*) for the point of the lowest pressure in the TMD line, p m 
is the point with the highest temperature and ph is the point with the highest pressure. 

The three top graphs in Fig.[2]show that the effect of decreasing the slope between the two 
length scales in the pair interaction potential is to move the TMD to higher temperatures. 
This result explains why the TMD is not observed in the discontinuous square well (DSW) 
model 8 . As the slope increases the TMD pressure and temperature approach the amorphous 
region and the system becomes unstable. For slopes higher than S3 case no anomalous 
behavior is observed. 

The middle graphs in Fig. [2] show that as the attractive scale increases, the TMD moves 
to higher temperatures and lower pressures as observed in potentials in which the attractive 
scale becomes dominant^. 

The bottom graphs in Fig. [2] show that as the repulsive scale becomes broader, the 
density anomaly region in the pressure temperature phase diagram goes to lower pressures 
and shrinks as observed in potentials in which the repulsive scale becomes dominant^. 

In addition in all the phase diagrams it is possible to observe that the TMD maximum 

9 



pressure never exceeds the critical pressure^. Our results indicate that the location in 
the pressure temperature phase diagram of the density anomalous region depends on the 
distance between the two length scales^"—. 
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FIG. 4. Diffusion coefficient as a function of density. The dots are the simulation data and the 
solid lines are polynomial fits. The dashed lines connect the densities of minimal and maximal 
diffusivity that limit the diffusion anomalous region. 



The Fig. H] shows the graphs of the dimensionless translational diffusion coefficient as 
function of density for all families, S, A and R. The solid gray lines are a polynomial fits 
to the data obtained by the simulations (the dots in the Fig. [3J. The diffusion coefficient 
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follows the same trend of the TMD line. This result is not surprising since the hierarchy of 
the anomalies suggests that the mechanism for the presence of a TMD line is related to the 
mechanism for the existence of a maximum and minimum diffusion coefficient. 

We also test the effects that changes in the two length scales have in the location in the 
pressure-temperature phase diagram of the structural anomalous region. 

The translational order parameter is defined as^ 1 ^ 1 ^ 



cutoff distance set to half of the simulation box times^ <?(£) is the radial distribution 
function which is proportional to the probability of finding a particle at a distance £ from a 
referent particle. The translational order parameter measures how structured is the system. 
For an ideal gas it is g — 1 and t — 0, while for the crystal phase it is g ^ 1 over long 
distances resulting in a large t. Therefore for normal fluids t increases with the increase of 
the density. 

The graphs in Fig. [5] illustrate the translational order parameter versus density for the 
potentials studied. The dot-dashed lines show the maximum and minimum in the values of 
t that limit the region of anomalous behavior. These extrema are also shown as dot-dashed 
lines in Figs. |2j The values at the pressure-temperature phase diagram for the different 
potentials follow the same trend as the TMD and diffusion anomalous regions. 

C. Radial distribution function 

The density anomaly can be related to the structure by analyzing the behavior of the 
radial distribution function. For a two length scales potential the g{r) has two peaks: one 
at the closest scale, r 1; and another at the furthest scale, r^. 

Recently it has been suggested that a signature of the presence of TMD line would be 
given by the radial distribution function as follows. At fixed temperature as the density is 
increased the radial distribution function of the closest scale, g(ri), would increase its value 
while the radial distribution function of the furthest scale, g(r 2 ), would decrease. This can 
be represented by the rule^ 1 ^ 




(4) 



where £ = rpl is the distance r in units of the mean interparticle separation p 3, £ c is the 



n lj2 



dgjr) 
dp 



X 



dgjr) 
dp 



< 



(5) 
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FIG. 5. The translational order parameter a function of density for fixed temperatures T* = 
0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5 (from top to bottom) for the families S (on top), A (on 
middle) and R (on bottom). The dot-dashed lines locate the maximal and minimal in t* . 

The physical picture behind this condition^ is that for a fixed temperature as density 
increases particles that are located at the attractive scale, r2, move to the repulsive scale, r%. 
Figures El illustrate a typical radial distribution functions for fixed T* as p* is varied. These 
graphs show that the picture of particles changing length scales due to pressure increase is 
valid for densities beyond a threshold density p* min - 

The regions identified by the radial distribution function as fulfilling the condition Eq. [5] 
are illustrated as opened circles in Fig. [7J The solid curve shows the TMD line. All the 

12 



stable state points with density equal or higher the minimum density at the TMD line verify 
the relation ITi )2 (p,T) < 0. This result gives support to our assumption that the presence 
of anomalies is related to particles moving from the furthest scale, r 2 , to closest length scale, 

n. 

The Figs. [8] show the value of the radial distribution function at the closest, g{r\) (dashed 
lines), and at the furthest scale, g(r 2 ) (solid lines), as a function of the reduced density, p* . 
For the closest scale g{r\) is monotonic with density while the value for the g(r 2 ) for a fixed 
temperature increases with the density for densities below the p* < p* min and decreases for 
densities above this threshold. This behavior, also shown in the Fig. [6j corroborates the 
condition stated in the Eq. [5] and supports the idea that particles move from one scale to 
the other by compression at p > p m in—- 

Besides to the move of particles from the attractive scale to the repulsive scale for p* > 




r r 



FIG. 6. Radial distribution as a function of the reduced distance for selected cases in the three 
families of potentials for T* = 0.8. For all the families, for p* < 0.40 the first and second peaks of 
g(r) increase with the increase of density. For p* > 0.40 the first peak increases while second peak 
decreases with the increase of density. 
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FIG. 7. Pressure-temperature phase diagram for the S, A and R families of potentials, illustrating 
as opened circles the regions where the condition Eq. [5] is obeyed. 

p* min as the pressure (density) is increased, particles also move from one scale to the other 
due to the increase of temperature^ for p* < p* min . At constant density, p* < p* min) the 
radial distribution function of the attractive scale, gfa), decreases with the increase of the 
temperature while g{r\) increases with the increase of temperature, indicating that particles 
move from one scale to the other due to thermal effects. At the density p* min1 the value of 
g{r\) is independent of the temperature. 

What is the meaning of the density p min in which g(ri) is independent of temperature? 
As it was pointed in the previous paragraph, for p* < p* min particles move from the furthest 
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p p p 



FIG. 8. g(ri) and gfa) for the S, A and R families of potentials as a function of the reduced 
density. The temperatures are T* = 0.4 — > 1.5 (from top to bottom p* > 0.6). The solid line 
connects the points, for different temperatures, where g{r\) = g{r2)- 



scale, r 2 , to the closest scale, r 1; using thermal energy as the temperature is increased. In 
this case g{r\) increases with temperature. For p* > p* min particles move from r 2 to 7*1, using 
the increase of pressure pressure (or density) as illustrated by the Eq. [51 From statistical 
point of view, the two mechanisms governing the behavior for p* > p* min and p* < p* min are 
quite different. While increasing temperature affects particles individually, increasing the 
density or the pressure affects the particles as clusters or networks. Then, as the potential 
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FIG. 9. Isobaric specific heat for the Si potential for different temperatures for different pressures. 

becomes more soft, the threshold density p* min beyond which the particles move form one 
scale to the other by compression should decrease as observed in Fig. El Therefore p* min is 
the threshold between these two mechanisms present in systems that have density anomaly. 

In addition to these low density limit, the density anomalous systems also have a high 
density threshold, p* max {T). Fig. E] illustrates as a solid thick line the temperatures and 
densities, Pmax(T), in which g(ri) = g{r 2 ). Since g(r) is related with the number of particles 
at distance r, for p* < p* max (T) more particles are in the attractive scale, r 2 , while for 
p* > p* max {T) more particles are at the repulsive scale r x . Therefore, the thick solid line 
is a boundary between the high density liquid and the low density liquid. Analysis of the 
stability indicates that no real phase transition is observed across this line. 

In order to understand what what happens in the region of the pressure-temperature 
phase diagram of the p* max {T), the behavior of the specific heat in this region was analyzed. 
Fig. [9] shows the curves of isobaric specific heat, for different pressures as a function of the 
temperature for the potential S%. The peak of c p for each one of the potentials analyzed in 
this manuscript coincide with the region p^^iT) where g{r\) = <?(r 2 ). This result indicates 
that the structure in the TMD region already build the liquid arrangements required for the 
liquid-liquid phase separation. 
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IV. CONCLUSIONS 



In this paper we have studied three families of core-softened potentials that exhibit two 
length scales, one repulsive, r 1; and another attractive, r 2 . 

We had observed that the region in the pressure-temperature phase diagram occupied 
by the TMD is quite sensitive to the slope between the two length scales. As the slope 
increases the region decreases. We also found that the region in the pressure-temperature 
phase diagram where the density, diffusion and structural anomalous behavior is observed 
shifts to lower pressures and shrinks the attractive scales or the repulsive scales become 
wider. Our results suggests the competition between two length scales are the relevant 
mechanism for the existence of the TMD. 

In an attempt to confirm this assertion we showed in this family of potentials that the 
condition 11^2 < seems to be associated with the presence of anomalous behavior. In addi- 
tion we also observed that the peaks of the radial distribution function at each length scales 
exhibit very distinct behavior with density and temperature suggesting two complementary 
mechanisms for the competition between the two scales. 

At low densities, p* < p* min particles move from r 2 to r\ with the increase of temperature, 
using thermal energy. For densities above p* min the increase of g(r{) is associated with the 
increase in the pressure (density). In this interval of densities at a certain density and 
temperature the radial distribution function at the first length scale equals the value at the 
second scale, namely g(ri) = g(r 2 ). This point can be identified with the Widom line. 

The relation between the radial distribution function and the Widom line, believed to be 
the onset of the liquid-liquid phase transition give the support to the idea that the Widom 
line separates two structural distinct regions that are also separated by a fragile-strong 
transition^. 

We expect that his result will not only shade some light in the definition of what is the 
shape an effective core-softened would have in order to held anomalies but also would serve 
to reinforce the ideas of linking dynamic transitions and thermodynamic properties. 
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